﻿#pragma once
#include <iostream>
#include <vector>
#include <eigen3/Eigen/Dense>
#include <fstream>
#include <cmath>
#include "Function.h"

using namespace std;
using namespace Eigen;

class Bform
{
private:
    vector<double> knots;
    vector<double> values;
    int N;
    int n;
    int condition;
    vector<double> result;

public:
    Bform(vector<double> _knots, vector<double> _values, int _n)
    {
        N = _knots.size();
        n = _n;
        double h = 1.0 * (_knots[N - 1] - _knots[0]) / (N - 1);
        for (int i = n; i > 0; i--)
        {
            knots.push_back(_knots[0] - i * h);
        }
        knots.insert(knots.end(), _knots.begin(), _knots.end());
        for (int i = 1; i <= n; i++)
        {
            knots.push_back(_knots[N - 1] + i * h);
        }
        values = _values;
    }

    Bform(Function& _A, vector<double> _knots, vector<double> _values, int _n, int _condition)
    {
        condition = _condition;
        N = _knots.size();
        n = _n;
        double h = 1.0 * (_knots[N - 1] - _knots[0]) / (N - 1);
        for (int i = n; i > 0; i--)
        {
            knots.push_back(_knots[0] - i * h);
        }
        knots.insert(knots.end(), _knots.begin(), _knots.end());
        for (int i = 1; i <= n; i++)
        {
            knots.push_back(_knots[N - 1] + i * h);
        }
        values = _values;

        if (condition == 0) {
            double m_a = _A.diff(_knots[0]);
            double m_b = _A.diff(_knots[N - 1]);
            complete(m_a, m_b);
        }
        else if (condition == 1) {
            double m_a = _A.diff2(_knots[0]);
            double m_b = _A.diff2(_knots[N - 1]);
            D2(m_a, m_b);
        }
        else if (condition == 2) {
            double m_a = 0;
            double m_b = 0;
            D2(m_a, m_b);
        }



    }
    double generate(int i, double _x, int _n)
    {
        if (_n == 0)
        {
            if (_x > knots[i - 1 + n - 1] && _x <= knots[i + _n + n - 1])
                return 1;
            else
                return 0;
        }
        else
        {
            return (_x - knots[i - 1 + n - 1]) / (knots[i + _n - 1 + n - 1] - knots[i - 1 + n - 1]) * generate(i, _x, _n - 1) + (knots[i + _n + n - 1] - _x) / (knots[i + _n + n - 1] - knots[i + n - 1]) * generate(i + 1, _x, _n - 1);
        }
    }

    double d(int i, double _x, int _n)
    {
        if (_n == 0)
        {
            return 0;
        }
        else
        {
            return (_x - knots[i - 1 + n - 1]) / (knots[i + _n - 1 + n - 1] - knots[i - 1 + n - 1]) * d(i, _x, _n - 1) + (knots[i + _n + n - 1] - _x) / (knots[i + _n + n - 1] - knots[i + n - 1]) * d(i + 1, _x, _n - 1) + 1.0 / (knots[i + _n - 1 + n - 1] - knots[i - 1 + n - 1]) * generate(i, _x, _n - 1) - 1.0 / (knots[i + _n + n - 1] - knots[i + n - 1]) * generate(i + 1, _x, _n - 1);
        }
    }

    double dd(int i, double _x, int _n)
    {
        if (_n == 0)
        {
            return 0;
        }
        else
        {
            return (_x - knots[i - 1 + n - 1]) / (knots[i + _n - 1 + n - 1] - knots[i - 1 + n - 1]) * dd(i, _x, _n - 1) + (knots[i + _n + n - 1] - _x) / (knots[i + _n + n - 1] - knots[i + n - 1]) * dd(i + 1, _x, _n - 1) + 2.0 / (knots[i + _n - 1 + n - 1] - knots[i - 1 + n - 1]) * d(i, _x, _n - 1) - 2.0 / (knots[i + _n + n - 1] - knots[i + n - 1]) * d(i + 1, _x, _n - 1);
        }
    }


    void complete(double m_0, double m_n)
    {
        MatrixXd A = MatrixXd::Zero(N + 2, N + 2);
        for (int i = 0; i < N; i++)
        {
            A(i, i) = generate(i - 1, knots[n + i], 3);
            A(i, i + 1) = generate(i, knots[n + i], 3);
            A(i, i + 2) = generate(i + 1, knots[n + i], 3);
        }
        A(N, 0) = d(-1, knots[n], 3);
        A(N, 1) = d(0, knots[n], 3);
        A(N, 2) = d(1, knots[n], 3);
        A(N + 1, N - 1) = d(N - 2, knots[n + N - 1], 3);
        A(N + 1, N) = d(N - 1, knots[n + N - 1], 3);
        A(N + 1, N + 1) = d(N, knots[n + N - 1], 3);
        VectorXd c(N + 2);
        for (int i = 0; i < N; i++)
        {
            c(i) = values[i];
        }
        c(N) = m_0;
        c(N + 1) = m_n;
        VectorXd R = A.lu().solve(c);
        for (int i = 0; i < N + 2; i++)
        {
            result.push_back(c(i));
        }
    }

    void D2(double m_0, double m_n)
    {
        MatrixXd A = MatrixXd::Zero(N + 2, N + 2);
        for (int i = 0; i < N; i++)
        {
            A(i, i) = generate(i - 1, knots[n + i], 3);
            A(i, i + 1) = generate(i, knots[n + i], 3);
            A(i, i + 2) = generate(i + 1, knots[n + i], 3);
        }
        A(N, 0) = dd(-1, knots[n], 3);
        A(N, 1) = dd(0, knots[n], 3);
        A(N, 2) = dd(1, knots[n], 3);
        A(N + 1, N - 1) = dd(N - 2, knots[n + N - 1], 3);
        A(N + 1, N) = dd(N - 1, knots[n + N - 1], 3);
        A(N + 1, N + 1) = dd(N, knots[n + N - 1], 3);
        VectorXd c(N + 2);
        for (int i = 0; i < N; i++)
        {
            c(i) = values[i];
        }
        c(N) = m_0;
        c(N + 1) = m_n;
        VectorXd R = A.lu().solve(c);
        for (int i = 0; i < N + 2; i++)
        {
            result.push_back(R(i));
        }
    }

    void natural()
    {
        D2(0, 0);
    }

    double cubic(double _x)
    {
        double s = 0;
        for (int i = 0; i < N + 2; i++)
        {
            s = s + result[i] * generate(i - 1, _x, 3);
        }
        return s;
    }


    double quad(double _x)
    {
        vector<double> t;
        for (int i = 0; i < N - 1; i++)
        {
            t.push_back((knots[n + i] + knots[n + i + 1]) / 2);
        }
        for (int i = 0; i < N - 1; i++)
        {
            t.push_back(((knots[n + i] + knots[n + i + 1]) / 2));
        }
        MatrixXd A = MatrixXd::Zero(N + 1, N + 1);
        A(0, 0) = generate(0, knots[n], 2);
        A(0, 1) = generate(1, knots[n], 2);
        for (int i = 1; i < N; i++)
        {
            A(i, i - 1) = generate(i - 1, t[i - 1], 2);
            A(i, i) = generate(i, t[i - 1], 2);
            A(i, i + 1) = generate(i + 1, t[i - 1], 2);
        }
        A(N, N - 1) = generate(N - 1, knots[n + N - 1], 2);
        A(N, N) = generate(N, knots[n + N - 1], 2);
        VectorXd c(N + 1);
        for (int i = 0; i < N + 1; i++)
        {
            c(i) = values[i];
        }
        VectorXd R = A.lu().solve(c);
        double s = 0;
        for (int i = 0; i < N + 1; i++)
        {
            s = s + R(i) * generate(i, _x, 2);
        }
        return s;
    }

};
